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I. INTRODUCTION 



Properties of strongly correlated quantum matter are 
usually well described by the many-body ground state 
and by the first elementary excitation. Multiparticle ex- 
citations are often not important because they just con- 
stitute an incoherent background. Thus, the study of 
quantum phase transitions mainly relies on low-energy 
spectrum analysis, namely energy levels of the ground 
state and the first excited state'^. However, in various 
systems, bound states may arise and play a major role. 
One of the most famous example are Cooper pairs which 
lead to superconductivity but other electron systems also 
display such mechanism (see, for instance, Refs. 2-4). In 
one-dimensional magnetic systems, well-known examples 
are the dimerized and frustrated spin chain as well as the 
two-leg spin ladder which contains bound states made 
up of triplons. Interestingly, such modes have been exper- 
imentally observed in cupratc ladder materials" . In two 
dimensions, the frustrated Shastry- Sutherland model and 
its experimental realization SrCu2 (603)2 are also known 
to possess two-triplon bound states^"''^ Note that more 
complicated bound states may arise in topologically or- 
dered systems where anyons (semions) can bind to form 
bosons or fermions as discussed in Ref. 14. 

Recently, such binding effects have been studied in the 
two-dimensional XXZ modeb ' where elementary excita- 
tions are dressed magnons. The aim of the present paper 
is to analyze the spectrum of such magnon bound states 
in two different spin systems. The first one is the fer- 
romagnetic transverse-field Ising model (TFIM) on the 
square lattice for which we derive the two-magnon spec- 
trum perturbativcly in the small- field limit. The high- 
order series expansion (order 12) of the corresponding 
gap allows us to extrapolate its behavior near the critical 
point where it is found to vanish. We also compute the 
ratio between the one-magnon gap and this two-magnon 



gap at the critical point which is approximately 1.8, in 
agreement with field-theoretical predictions^'"'". The 
second system considered in this study is the XXZ model 
on the square lattice which, in the isotropic limit (XXX), 
is the celebrated antiferromagnetic Heisenberg model. As 
for the TFIM, we focus on the two-magnon bound-states 
spectrum which is derived up to order 8 near the Ising 
limit. However, as we shall see, results obtained here 
differ from those obtained recently by Hamcr' ' since, 
contrary to his claim, we show using very simple ar- 
guments, that there are two distinct branches of two- 
magnon bound states at low energy. Also, let us under- 
line that the fate of the lowest-energy bound state when 
approaching the Heisenberg limit cannot be conclusively 
determined at this order. 

From a methodological point of view, several meth- 
ods to compute high-order series expansion in quan- 
tum many-body systems are available (see Refs. 18,19 
for a review). Here, we use the perturbative con- 
tinuous unitary transformations (PCUTs) method-*^"^'^ 
which is especially well suited to investigate the many- 
particle spectrum and provides a natural quasiparticle 
(QP) description. 

The structure of this paper is the following. In Sec. II, 
wc introduce the two models (TFIM in Sec. II A and XXZ 
model in Sec. II B) under consideration. Wc show that 
the Ising limit is a good starting point for a perturbation 
theory in the PCUTs framework and we also give a very 
simple picture to understand the occurrence of bound 
states in this limiting case. In Sec. II C, we introduce an- 
other description of the spin model in terms of bond de- 
grees of freedom which will be useful to set up the present 
perturbation theory framework. The end of this Sec. II D 
is dedicated to symmetry considerations. In Sec. HI, we 
recall several important aspects of PCUTs which are es- 
sential to understand the next section. In Sec. IV, we 
adapt the finite-lattice method (commonly used in statis- 
tical mechanics "^'^^) to quantum problems, allowing one 
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FIG. 1: (Color online) Square lattice on which the models 
studied in this paper are defined. Translation vectors in the x 
and y directions are denoted n\ and n^. The bond between 
sites i and j is denoted or by a Greek letter such as /3 

to shorten the notation. 



FIG. 2: (Color online) Two magnons (spins pointing left) 
on top of a ferromagnetic background. Such a state involves 
eight antiferromagnetic (bold) bonds (four for each magnon), 
and costs an energy 8, with respect to Hamiltonian Hi. 



to significantly increase the maximum order of the series 
expansions. Let us stress that readers interested only in 
results can skip these two sections and switch directly to 
Sees. V and VI where we discuss the low-energy spectra 
of the TFIM and XXZ model, respectively. Finally, in 
Sec. VII, we discuss the spectrum of a new model which 
naturally emerges when describing the TFIM in terms 
of bonds. This model may be seen as a special case of 
the toric code model" ' in a magnetic field'^ in which 
flux creation energy cost vanishes. All coefficients used to 
compute gaps in both models are gathered in Appendices 
A and B. 



II. MODELS AND MAPPINGS 

A. Transverse-field Ising model 

Let us first introduce the TFIM on the square lattice 
whose Hamiltonian reads 

i^TFIM ^-jY^ <^>] ~hY<^'l. (1) 

where the second sum runs over all i of the square lat- 
tice and the first sum runs over all bonds between 
nearest- neighbor sites (see Fig. 1). The Pauli matrices 
are denoted by ti^, cr^, and cr^. Here, we focus on a 
ferromagnetic Ising coupling J > and, without loss of 
generality, we consider ft, > 0. 

This model whose classical counterpart is the three- 
dimensional Ising model is known to display a second- 
order phase transition separating a symmetric (disor- 
dered) phase at large field h from a broken (ordered) 
phase at large coupling J. For an infinite field, all spins 
point in the -Vx direction, while for an infinite exchange 
coupling, all spins point either in the -\z direction or in 
the —z direction. 



The ground-state energy as well as the single- 
excitation spectrum have been computed in both phases 
using series expansion^"--' '-'^ allowing for a precise de- 
termination of the critical point at J//i|c = 0.3285(1). 
In the broken phase, the perturbative expansion is done 
around the Ising limit (h ~ 0). Thus, it is convenient to 
introduce the following Hamiltonian that will also be the 
starting point of our study of the Heisenberg model (see 
Sec. HB) 

ii. = -\Y.-lo]. (2) 



Let us consider its ferromagnetic ground state where 
all spins point in the -\-z direction. Lowest-energy exci- 
tations then consist in static magnons (spin flips) whose 
energy cost is 4. The prefactor 1/2 in Eq. (2) defining H\ 
indeed ensures that any state will have an energy equal to 
that of the ground state plus the number of antiferromag- 
netic bonds. Two such "isolated" magnons are shown in 
Fig. 2. 

A two-magnon state will have an energy cost of 8, ex- 
cept if the two magnons are nearest neighbors. Indeed, 
in such a case, only six bonds have an antiferromagnetic 
configuration, which results in an energy reduction of 
8 — 6 = 2 compared to the situation where magnons are 
far apart. This proves the existence of bound states in 
the spectrum, two of which are shown in Fig. 3. 

Of course, following the same line of reasoning, one can 
show that n-magnon bound states exist for any n ^ 2, 
but we shall restrict our study to the case n = 2. Let us 
also notice that in the ordered phase, the perturbation is 
performed around the field term of Eq. (1), which basi- 
cally counts spin flips, so that no binding effect is present 
in this limit. 
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FIG. 3: (Color online) Two examples of two-magnon bound 
states. Each of them has an energy cost of 6, with respect to 
Hamiltonian Hi. As in Fig. 2, antiferromagnetic bonds are 
represented as bold segments. 



B. Heisenberg and XXZ models 

Let us now discuss, in the same spirit, the antiferro- 
magnetic Heisenberg model on the square lattice. Its 
Hamiltonian reads 



(3) 



where Si — cri/2 is the spin operator at site i. In the 
following, we set J = 2 and we introduce an anisotropy 
parameter A with the aim of performing series expansion 
in this parameter, as was done in Refs. 2!) and 1!),;^0,-']1. 
The Heisenberg Hamiltonian is then the A = 1 limit of 
the following XXZ Hamiltonian 



where wc have introduced the usual raising and lower- 
ing operators = ^{cr'fzLiaf). The square lattice 
being bipartite, it is possible to perform a rotation of 
angle tt around the x axis on one of the sublattices. 



namely (cr'^,CT^,(T^ 



This transfor- 



mation leads, without changing notations for the Pauli 
operators and making use of Hi defined previously, to 



XXZ 



= Hi + X 



(6) 



Wc thus sec that in the limit of a vanishing anisotropy 
parameter A, Hxxz reduces to Hi so that the binding 
effects discussed above for -ffTFiM arc also at work here. 
The difference between both models is the perturbation 
term added to Hi. 

Let us finally mention that from the point of view of 
the XXZ model, a second-order quantum phase transi- 
tion occurs at the Heisenberg point A = 1, separating a 
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FIG. 4: (Color online) Dots are located in the middle of the 
bonds of the original square lattice of Fig. 1 (represented with 
dashed lines). Bond operators of Eq. (7) are defined on these 
bonds. White plaquettes (centered on vertices of the original 
lattice) are used to define Aa operators [see Eq. (9)] , whereas 
Bp operators of Eq. (12) are defined on gray plaquettes. The 
bond that is shared by two white plaquettes such as s and t is 
denoted by (s, t) or by a Greek letter such as a. Contours Ci 
and C2 are used in Eq. (14) to define Z2 conserved quantities 
for a system with periodic boundary conditions. 



gapped phase for < A < 1 with a twofold-degenerate 
ground state with Neel order from a gapless phase with 
0(2) symmetry for A > 1. 



C. Bond description 

As we have already seen in Sec. HA when describ- 
ing the spectrum of Hi, the energy cost of any multi- 
magnon state is given by the number of antiferromagnetic 
bonds of the state's spin configuration. It is thus natu- 
ral, though not mandatory, to introduce effective spin 
variables living on the bonds (see Fig. 4) as follows : 



(7) 



To make notations light, wc denote bonds between 
two sites i and j with Greek letters, such as /? in the 
above equation. A value -1-1 or —1 of ct^ is then asso- 
ciated to, respectively, a ferromagnetic or antiferromag- 
netic configuration of bond /3. Within such a description, 
Hamiltonian (2) becomes a pure field term 



Hi 



(8) 



/3 



We can now give another form of the field term of 
ffxFiM- Since erf flips the spin at site i, it flips the four 
bonds sharing site i. Let us denote by s{i) the set of these 
four bonds (the notation s referring to stars or vertices 
of the original lattice). Furthermore, we introduce 



A, 



n 



(9) 
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which is the product of the four operators flipping bonds 
connected to site i. One can alternatively say that the 
As are defined on the white plaquettes of Fig. 4. Then, 
the TFIM Hamiltonian reads 

i?TFIM = -i$]a|-/i^I.. (10) 

In the same vein, we can rewrite the XX part of the 
XXZ Hamiltonian (4). From expression (6), it is clear 
that this XX term flips two adjacent spins, if these two 
spins are in a ferromagnetic configuration while it annihi- 
lates antiferromagnetic bonds. We recall that in Eq. (6), 
the sublattice rotation was already used, so that a fer- 
romagnetic configuration corresponds to an original an- 
tiferromagnetic configuration in Eqs. (4) and (5). As a 
consequence, the XX term flips the six bonds around a 
ferromagnetic bond, and one can write 

i^xxz = -IJ2^i + ^J2 AsM^f^, (11) 

/3 {s,t) 

where the second sum runs over nearest-neighbor white 
plaquettes s and t that share a common bond a = (s,^) 
(see Fig. 4), and involves projectors that annihilate an 
antiferromagnetic configuration of such bonds. 

D. Symmetries, counting, and relation to toric 
code 

Up to now, we have not yet discussed one fundamen- 
tal aspect of the bond description, namely, the counting 
of states. Indeed, if we assume that the original lattice 
contains N spins half cr, the mapping brings us to a de- 
scription with 2N spins half 5" since the number of bonds 
is twice the number of sites (assuming periodic boundary 
conditions). Thus, it seems that Hilbert spaces of the H 
Hamiltonians (8), (10), and (11) are too large, and in- 
volve "unphysical" states. 

However, looking at these Hamiltonians a bit closer, it 
is clear that the following operators are conserved 

Bp=\{a^p. (12) 

In this definition, the product is performed over all bonds 
belonging to p, which can be any of the gray plaquettes 
shown in Fig. 4 (they are also plaquettes in the original 
lattice shown in Fig. 1). Among the TV I2 operators that 
can be defined this way, only iV — 1 can be set indepen- 
dently to ±1, because of the constraint J^^ Bp ^ 1 (with 
periodic boundary conditions). 

With such notations, Hamiltonian (10) is nothing but 
the toric code Hamiltonian- ' in a magnetic field (we keep 
using the ~ notations) 

s p p 



with Js = /i; Jp = and hz = J = 1/2. Although 
the original toric code (with nonvanishing Jg and Jp) in 
a field has already been the subject of some works (see 
Refs. 32, and 26, .33, 34), we are not aware of any study of 
this precise model which has also been recently obtained 
in a related framework ' '. For this reason, we will discuss 
some of its properties in Sec. VII. 

For periodic boundary conditions, as in the "bare" 
toric code"'"' (see also Ref. 3() for a pedagogical descrip- 
tion), -Bp's are not the only conserved operators when a 
magnetic field in the z direction is present. Using con- 
tours Ci and C2 depicted in Fig. 4, one can indeed define 
the following operators 

Ci = [] and C2 - n ^1- (14) 

/3eCi /3eC2 

These are also Z2 conserved quantities, which can be set 
to ±1 independently of the values of the Bp's. 

All in all, we have {N — l)-|-2 = iV + 1 conserved 
and independent Z2 quantities. Furthermore, to recover 
the "physical" subspace and the physics of the TFIM or 
the Heisenberg model, one should set all of them to -t-1. 
This reduces the Hilbert space dimension from 2^^^ to 
22W-(Ar+i) ^ 2^-1, that is in fact less than the original 
Hilbert space dimension. However, the reason for this ob- 
vious : our description in terms of bond variables is done 
with respect to one of the two degenerate ferromagnetic 
ground states. As a consequence, such a description can 
only be valid in the broken phase. 

III. PERTURBATIVE CONTINUOUS UNITARY 
TRANSFORMATIONS 

The aim of this section is to equip the reader with the 
basic knowledge about PCUTs necessary for the under- 
standing of the next section. For concreteness, we work 
out an example in detail, rather than focusing on a gen- 
eral framework. 



A. Basic ideas of continuous unitary 
transformations 

The CUTs method as known in the condensed- 
matter theory community originates from the work of 
Wegner^"'"*' A pedagogical introduction to this tech- 
nique can be found in Refs. 3!) and 40. The aim of this 
technique is to diagonalize or, more modestly, to block 
diagonalize a given Hamiltonian H thanks to a unitary 
transformation. The latter is not performed in a single 
step but rather in a continuous way (whence the name of 
the method) as 

H{l) = U\l)HU{l), (15) 

where I is a running parameter such that H = H{1 — 0) 
and iJoff = H{1 = 00) is an effective (block-) diagonal 



5 



Hamiltonian. This equation can be cast into a differential 
(flow) commutator equation^" 

diHil)^m,Hil)], (16) 

where r]{l) — diW {l)U{l) is the anti-Hermitian generator 
associated to the unitary transformation U{1). 

B. Quasi-particle conserving generator 

The next task is to find the appropriate generator 
which, from the local knowledge of H{1), leads to the 
desired form of the effective Hamiltonian. We shall only 
discuss the QP conserving generator-- that will be used 
in the sequel. We furthermore focus on a specific example 
for which the Hamiltonian can be written 



H = 



E 



Tn 



(17) 



For concreteness, in the following we consider the case 
where n £ {0, ±2, ±4} which is relevant for the TFIM. In 
this equation, Q is the Hermitian operator which counts 
the number of QPs (so its spectrum is contained in N), 
and the T„ 's are operators that change the QP number by 
the amount n, so that [QjTn] = 7iT„. The hermiticity of 
the Hamiltonian requires that T,t = T_„. The QP con- 
serving generator is designed to bring the Hamiltonian 
to an effective form that conserves the number of QPs : 
[Q, Hcs] = 0. Said differently, under the CUTs, aU terms 
Tn will be fiowing (but not Q which is isolated from all 
other terms), thus becoming Tn{l), and one wishes to 
reach a situation where T„(Z = oo) = for all n ^ 0. 
This can be achieved--' by choosing 



77(0 = r+2(0 - T_2(0 + r+4(0 - t_4(0- 



(18) 



With this choice of generator, the flow Eq. (16) can be 
written 

9,To(0 = 2[r+2(0,r-2(0] +2[T+4(/),T_4(0], 
diT+2{l) = -2T+2{l)+[T+2{l),n{l)\ 

+2[r+4(0,r_2(0], 
diT+S) - -4r+4(0 + [r+4 (0,7^0(0]- (i9) 

We have not written the flow equations for T_2(/) and 
T-i{l) since the Hamiltonian remains Hermitian under a 
unitary transformation. Let us emphasize that no new 
term appears during the flow, thanks to the choice of 
the generator (e. g., no term creating six particles ap- 
pears since [T+2(/),T+4(/)] + [T+4{l),T+2{l)] = 0). Lin- 
ear terms in the right-hand side of Eq. (19) ensure that 
T„^o(? = cx)) = so that [iJcff , Q] = 0. 



C. Perturbative commutator expansion of the flow 
equation 

It remains to solve these flow equations which is still a 
challenging task since the Tn{l) terms contain an infinite 



number of operators. The easiest way of doing so is to 
perform a perturbative expansion of the flow equations, 
assuming that all r„ terms in H [see Eq. (17)] are small, 
and of the same order of magnitude. Such an expansion 
was first performed by Stein- ^ for a Hamiltonian with- 
out r±4 terms, and extended to Hamiltonians with terms 
changing Q by any amount by Knetter and Uhrig-- who 
also provided a description in terms of QP. However the 
formalism used in these papers was not the same as the 
one presented here, and, in particular, the emphasis was 
not on getting a commutator expansion, which is our goal 
in the following. 

To this end, the expansion for T„(Z) is simply written 



r«(0 = E^»HO, 



(20) 



where the superscript (i) is the order, in perturbation, of 
Ti'\l). The flow equations can then be expanded as 



+2 [Tiqil),T[^r\l)]) , 

+2 [Tiqil),T^_^P\l) 



ftr«(0 



3 = 1 



These have to be solved with the initial conditions 
Tn\l = 0) = di^iTn and then one can take the limit 
Z — > oo to obtain iJoff- For the Hamiltonian considered 
up to now, one obtains, at order 3, the following commu- 
tator expansion 



^^cff — 



Q + ro + ^[r+2,T_2] + J[r+4,T_4] (22) 



32 



[[T+2,To],T_2] + [T+2,[ro,T_2]] 



[[T+4,T_2],T_2] + [T+2,[r+2,T_4]] 
[[r+4,To],T_4] + [T+4,[To,r_4]] 



We have computed such expansions for Hamiltonians 
of the form given in Eq. (17) to various orders given in 
Table I. In this table, we also give the total number 
of nonzero coefficients one obtains once the commuta- 
tors in the effective Hamiltonian have been expanded in 
polynomials in T operators. Let us note that for the 
Hamiltonian we have been considering up to now, with 
n G {—4, —2, 0, 2, 4}, the expansion can be obtained from 
the one of rimax = 2 (i e., n G {—2,-1,0,1,2}) by a 
proper rescaling of the corresponding coefficients. 
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Let us also mention that all coefficients are given by 
rational numbers and are valid for arbitrary system size, 
including the thermodynamical limit. The very last (but 
not least) step in the whole computation is to apply this 
effective Hamiltonian to states with fixed number of QPs 
and to diagonalize it. 

To determine the low-energy QP properties of the 
TFIM and the XXZ Hamiltonian around the Ising limit, 
we shall use this PCUTs approach. Indeed, in this limit, 
both systems can be written in the form (17) by defining 
the operator 



E 

/3 



N + Hi, 



(23) 



which counts the number of antifcrromagnetic bonds. 
The operators T„ arc then proportional to the perturba- 
tion {h in the TFIM and A in the XXZ model). The index 
n denotes the change in the number of antifcrromagnetic 
bonds q. For the TFIM one has n e {0, ±2, ±4} and 
for the XXZ model one gets n e {0, ±2, ±4, ±6}. The 
larger number of T„ operators for the XXZ model results 
in a larger effort since more processes have to be taken 
into account for a given perturbation order. For instance, 
for the 2QP sector (containing bound states), we reached 
order 12 for the TFIM and order 8 for the XXZ model. 

Here, we use the transformed Hamiltonians (10) and 
(11) defined with bond variables, Eq. (7), but PCUTs 
can be (and for actual computer implementations arc) 
applied directly to Eqs. (1) and (6). In order to avoid pos- 
sible confusion, we shall use the notation OQP, IQP, and 
2QP when referring to 0, 1, and 2 magnons, and Oqp, Iqp, 
. . . , 8qp when referring to 0, 1, ... ,8 antifcrromagnetic 
bonds. The ground state of the effective Hamiltonian is 
the state without antifcrromagnetic bonds (q = 0). The 
one-magnon excitations have q = A and bound states of 
two nearest-neighbor magnons have q = 6 antifcrromag- 
netic bonds. Things become much more complicated for 
configurations with more antifcrromagnetic bonds, and 
we shall not study them here. Let us simply mention 
that q = 8 corresponds, in the unperturbed limit, either 
to two unbound magnons or to three- or four-magnon 
bound states. However, at finite coupling, we cannot 
exclude that there also exist two-magnon bound states 
(different from those discussed above) in this sector. 





Order 


Number of coefficients 


1 


18 


67214380 


2 


14 


569842124 


3 


12 


924457284 


4 


10 


189956506 



TABLE I: Maximum order at which we have derived the effec- 
tive Hamiltonian, for various values of rimax [see Eq. (17) for 
definitions] , as well as the total number of non-zero coefficients 
needed to express the effective Hamiltonian as a polynomial 
in r„'s. 



IV. LINKED-CLUSTER EXPANSION AND 
QUANTUM FINITE-LATTICE METHOD 

A. Linked-cluster expansion 

Flow Eq. (21) show that iJoff can be written as a per- 
turbative commutator expansion to any order, as exem- 
plified in Eq. (22) to order 3. This property has dramatic 
consequences when one considers a Hamiltonian defined 
on a lattice with local r„ operators. By local, we mean 
that Tn = where i runs over the lattice sites, 

with Tn.i acting on a finite number of sites neighboring i. 
Indeed, in such a situation, commutators [T„^i,Tp van- 
ish as soon as i and j arc sufficiently far apart and are 
local operators. One could try to implement the calcula- 
tion of the commutators in a symbolic way but it is usu- 
ally much easier to apply the effective Hamiltonian (with 
expanded commutators) to states. For the practical pur- 
pose of evaluating the action of one term in iJoff, one then 
only needs to apply Haft to a finite-size linked cluster of 
sites. The notion of linked is problem-dependent since it 
depends on the extension of T„ ^ operators. 

The appearance of a linked-cluster expansion for ef- 
fective Hamiltonians written as commutators is a long- 
established and well-known result (see e. g., Ref. 41), al- 
though it seems to have escaped the attention it deserves 
in more recent publications'"^'^"'" '. 

Let us finally note that the PCUTs method is not the 
only available one to obtain such a pcrturbativc commu- 
tator expansion and one could use the Van Vlcck for- 
malism (sec e. g., Ref. 42) as well. The main advantage 
of PCUTs, is that it does not generate terms creating 
any number of QPs. so that the bookkeeping is not too 
difficult. One drawback is that one has to solve flow 
equations instead of performing purely algebraic manip- 
ulations. 

In the following, we will introduce the two com- 
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FIG. 5: (Color online) Pictorial representation of the ground- 
state energy per site as a linked-cluster expansion. A cross 
on a site means that this site has been acted upon at least 
twice by a T„ operator (two spin flips per site are necessary 
to preserve the ferromagnetic state). The second line follows 
from the symmetries of the lattice. All contributions relevant 
for a computation at order 3 of eo are shown, contrarily to 
contributions of order 4 (which are all gathered in ". . ."). 
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FIG. 6: (Color online) Pictorial representation of the compu- 
tation of matrix elements of HbH for various clusters. Crosses 
have the same meaning as in Fig. 5. Use of symmetry is made 
to simplify equations and to compute only four graphs (see 
Fig. 5) but others can be computed as easily. 



monly used implementations of the linked-cluster theo- 
rem which represent extreme cases. Either a calculation 
is performed on a large number of minimal clusters or a 
calculation is done on one very large cluster. Afterwards, 
a third implementation of the linked-cluster theorem is 
presented which is a compromise and which we use in the 
current study. 



B. Graph-ology 

With Sec. IV A in mind, it should now be clear that the 
aim is to compute contributions of the effective Hamilto- 
nian when the latter is allowed to act on states defined on 
clusters of one site, two sites, etc. To make things more 
concrete, in what follows, we shall focus on the computa- 
tion of the ground-state energy per site eo of the TFIM at 
order 6 in /i. As only even powers of h will appear in the 
expansion, we shall for simplicity only refer to the order 
in in this section, so that we will say that we work at 
order n when computing terms up to order (h?)^ . Only 
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FIG. 7: (Color online) Pictorial representation of the in- 
version of the relations of Fig. 6, which leads to subcluster 
subtraction. The last equality is due to symmetries. 



the action on the ferromagnetic (OQP) state will be con- 
sidered here so that we do not draw any arrows to make 
things clearer. The relevant contributions from different 
clusters are shown in Fig. 5. Let us mention that the 
clusters shown in this figure do have contributions at or- 
der 4 or more (in which case at least one site is acted 
upon twice) but clusters not shown (with four crosses or 
more) do not have contributions at order 3 or less. To 
obtain the last line of Fig. 5, the symmetries of the lattice 
have been used. Note that the r„ operators can act on 
at most three sites (each pictogram has at most 3 crosses 
in Fig. 5) at order 3 and in the OQP sector, because each 
r„ operator flips the spin of the site on which it acts, 
so one needs two r„ operators per site to start from and 
end up with a ferromagnetic (OQP) state. 

To compute the contributions shown in Fig. 5, it is 
not very practical to make sure that all sites have been 
acted upon at least once by a r„ operator. It is much 
easier in actual computations to compute all possible out- 
comes of the action of H^s on a given cluster and to sub- 
tract contributions of subclusters, as illustrated in Figs. 6 
and 7. Let us note that these subtractions are manda- 
tory in the perturbation theory used in Refs. 18 and 19 
which is not based on an expression of Hctt as a series in 
Tn operators. Numerical coefficients appearing in Fig. 7 
arc simply the number of ways the subclusters (or their 
symmetric-related ones) can be embedded in a given clus- 
ter. 

Thus, we are naturally led to enumerate all possible 
linked clusters [i. e., also graphs, hence the name of this 
section) that can be embedded in the square lattice, and 
apply the subtraction technique to each of the cluster. 
Though this is the least memory-consuming way to pro- 
ceed, and though the computational effort required for 
the computation on one cluster is rather small, it re- 
quires to perform a heavy and time-consuming combina- 
torial work. One indeed has to enumerate all possible 
linked clusters, and then, for each cluster, one needs to 
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FIG. 8: (Color online) Fourth (and last) matrix element from 
the 2x2 rectangular cluster to be considered for a finite-lattice 
computation at order 3 (apart from the first three of Fig. 6). 
The term in parentheses involves the action on four sites and 
only appears at order 4 and beyond (whence the exponent of 
the parenthesis). 



find all its linked subclusters. In order for this technique 
to be as efficient as possible, one should also ensure that 
topologically identical clusters are identified (such as the 
two three-site clusters of Fig. 7) . This again reduces the 
computational effort but makes the combinatorial tasks 
even harder and more time consuming. 

A completely opposite way of performing the calcula- 
tion is to reduce the combinatorial complexity to zero by 
computing all contributions in one go, thanks to a clus- 
ter with periodic boundary conditions, large enough so 
that it can accommodate all clusters one is interested in, 
without any finite-size effect. In this way, applying H^s 
to a OQP state defined on such a cluster one recovers the 
same state, up to a multiplicative factor being equal to 
Co times the number of sites of the cluster. However, as 
the Hilbert-space size quickly grows with the number of 
sites, this makes all computations greedy in memory, but 
also in time (because the number of intermediate states 
generated when applying Hcs can become huge). 



C. Quantum finite- lattice method 

We shall now see how one can work half-way between 
these two extreme cases. The idea is to generalize Ent- 
ing's finite-lattice method (see Refs. 19,24,4.3, and refer- 
ences therein) from the statistical physics setting (where 
it is used to compute the free energy) to the realm of 
quantum physics. Essentially, the idea is to use rectan- 
gular clusters only and to perform appropriate subtrac- 
tions. The main advantage of using rectangular clusters 
is that computing the number of cmbeddings of rectan- 
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FIG. 9: (Color online) Inversion of the equality of Fig. 8, 
with the help of the first three lines of Fig. 7, showing how 
subcluster subtraction appears. The second equality makes 
use of symmetries, contrarily to the first line. 



gular subclusters is trivial. Compared to graph-ology, 
one gains enormously from the combinatorial side, but of 
course one has to pay a price : the clusters one uses are 
larger than those of graph-ology. Enting's finite-lattice 
method is most useful in low dimensions (namely, two di- 
mensions in statistical mechanics) where the graph sizes 
do not grow too fast with the order in perturbation but 
with modifications it can also be efficiently applied to 
higher dimensions (sec Ref. 44 for an application in three 
dimensions). 

Let us first discuss this method in the OQP sector, 
where it is almost the same, in the context we arc work- 
ing in, as the original method of statistical mechanics, 
because the ground state is not infinitely degenerate. In 
the latter case, our method would provide an effective 
Hamiltonian in the low-energy subspace which is a major 
difference compared to a single number (the ground-state 
energy or the free energy in statistical mechanics). The 
first three equations of Fig. 6 remain valid since they in- 
volve a rectangular cluster but the last equation is now 
replaced by the new matrix element of Fig. 8. 

This can again be inverted to yield subcluster sub- 
traction and is shown in Fig. 9. From this, one under- 
stands that a rectangular cluster yields, after subtraction, 
the sum of all contributions of the linked subclusters that 
cannot be embedded in any rectangular subcluster ( of the 
considered cluster). If the cluster has width or height 1, 
the sum actually consists in a single term. From these 
considerations, one can deduce a simple formula. Let us 
call Ds^^s the expectation value of i/cff, computed on a 
rectangular cluster of size (s^, Sy), so that the first three 
equations of Fig. 6 and the one of Fig. 8 correspond to 
Di^i, £'2,1, £'3,1, and D2.2- Let us furthermore define 
(recursively) 



,,Sy ^Ds^^Sy-}^ iSx-t^ + l){Sy-ty + l)Dt^ 



(24) 



where means that the sum is restricted to the set 
of strict subclusters of sizes {tx,ty) of the cluster of size 
(sxiSy), namely, satisfying I ^ t^ ^ Sx and 1 ^ ty ^ Sy 
with (tx, ty) 7^ {sx, Sy). As an illustration, the first three 



lines of Figs. 7 and 9 give -Di,i, -D2,i, D^^i, and I?2,2- 
Then the ground-state energy can be expressed as a sum 
of Dg^^sy contributions, for a range of Sx and Sy that is 
problem dependent (because the T„ operators can have 
different spatial extensions, as is the case for the TFIM 
and the Heisenberg model). In the case of the TFIM we 
have been focusing on up to now, one has 

eo = e(°)+^e("\ (25) 

n 

with 

e("^i)= ^ 5,^,,^, (26) 

l^s^+Sy-liin 

and Cq^'' is the constant shift of order in the Hamil- 
tonian. Contrary to the TFIM where T„ operators act 
on one site, the r„'s for the Heisenberg model act on 

(n) 

two sites so that the sum for Cq would be restricted to 
I Sx + Sy ~ 1 ^ 2n, where n is the order in A^. 

The very same idea can be applied for any QP sector, 
provided appropriate subtractions of contributions from 
sectors with lower number of QPs are performed, as in 
standard linked-cluster expansions^'""'^'''^''. In the same 
spirit as before, let us give a pictorial explanation based 
on the TFIM, at order 2, and for the IQP sector (the 2QP 
sector can be worked out in a similar way) . The aim is to 
compute hopping amplitudes, as those given in Appendix 
Al. As can be seen in Fig. 10, one difference compared 
to the OQP sector is that depending on which term of 
the Hamiltonian one computes, the particle can hop to 
different final positions. In the figure, we have again not 
represented any arrow to make things light, but it should 
be clear that the reference state is a ferromagnetic state, 
with one spin flip at the QP's position. We have only 
considered hoppings tij of i sites in x direction and j 
sites in y direction, for which i ^ j ^ 0, other processes 
can be found thanks to symmetries of the square lattice 
and the hermiticity of the Hamiltonian. The hopping 
amplitude io.Oj which should rather be called a chemical 
potential, is a bit peculiar. First, the counting operator 
Q appearing in HcS gives a contribution of 4. Second, as 
one is only interested in the excitation energy above the 
ground state, the OQP contributions of the clusters have 
to be subtracted. 

Various contributions of Fig. 10 can be extracted from 
the effective Hamiltonian's matrix elements which are 
shown in Figs. 11 and 12. Again, one can invert these 
relations, which leads to a recursive subcluster subtrac- 
tion scheme (see Figs. 13 and 14). It is not as easy as in 
the OQP sector to write down a concise formula similar 
to Eqs.(24) and (26), because the particle position has 
to be taken in consideration, and implies some more re- 
strictions. However things should be clear from Figs. 11 
and 12. In practice, it is much easier to use a computer 
program that determines all relevant subclusters and sub- 
tracts their contributions. 
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FIG. 10: (Color online) Contributions to hopping amplitude 
at order 2, as a linked-cluster expansion. Crosses have the 
same meaning as in Fig. 5, and are acted upon at least twice 
by T„ operators. The empty (full) circle represents the initial 
(final) position of the particle. Sites with circles are acted 
upon at least once by T„ operators. When the "particle" does 
not move, i. e., for to,o, the particle's site (only represented 
with a full circle) is acted upon at least twice by T„ operators. 
In this case, as one is only interested in the energy of a IQP 
state with respect to the OQP ground state, one must subtract 
the OQP amplitude (Refs. IS, 19 and 2.3). Note that for to,o, 
one should also not forget the action of the counting operator 
Q, which gives a zeroth-order contribution equal to 4. 

As a final word, let us mention that the 2QP sector 
can be treated in the same way since the bound states in 
the TFIM or the Heisenberg model, that are of interest 
to us, behave as a single particle (they are made of two 
nearest-neighbor magnons). If one was interested in the 
2QP sector where the two magnons can move freely, one 
would not only need to perform the OQP subtraction as 
previously, but one would also have to subtract the two 
IQP contributions, in order to extract the 2QP interac- 
tions. We shall now apply this formalism to both systems 
introduced in Sec. II. 



V. RESULTS FOR TFIM 

A. Series expansion of the low-energy spectrum 

As we have seen, PCUTs provide, order by order, an ef- 
fective Hamiltonian Hj^^^^ which is unitarily equivalent 
to i?TFiM and commutes with the operator Q counting 
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FIG. 12: (Color online) Matrix elements of Hcn between 
different initial and final states that are useful for the com- 
putation of tifi, t2,o, and ti.i at order 2. Conventions are the 
same as in Fig. 11. 



FIG. 11: (Color online) Matrix elements of HcH between 
identical initial and final states that are useful for the compu- 
tation of tofi at order 2. Diamonds with solid lines represent 
the particle position in the considered state. Contributions 
arising from an action of H^s — Q which does not involve the 
particle site contain a dotted diamond that keeps trace of this 
site. Other symbols have the same meaning as in previous fig- 
ures. Terms in parenthesis are higher-order terms, for which 
the exponent gives the lowest order they start to contribute. 
For the last contribution, intermediate steps have not been 
shown, and the higher-order terms are not kept until the end, 
to keep things as light as possible. 



the number of antifcrromagnetic bonds. One must then 
analyze Hj^^^^ in each sector with a given number q of 
antiferromagnetic bonds. In the present work, we de- 
rived this effective Hamiltonian in the sector q = (Oqp 
or OQP) up to order 14 while we reached order 12 for 
(7 = 4 (4qp corresponding to IQP) and q = 6 (6qp corre- 
sponding to a bound state of 2QP). We restrict ourselves 
to the "physical" subspacc {i.e., the states correspond- 
ing to magnons) where all conserved quantities (i?p's) of 
Sec. II D are set to -1-1, and postpone the discussion of 



"unphysical" states (having antiferromagnetic bonds not 
corresponding to a magnon configuration of the initial 
models) to Sec. VII. 

For q ~ 0, wc obtain the ground-state energy per bond 



6144 



en = h h 

" 2 8 384 

1388129 10 

"254803968000 

707258321166713 

"2588971389419520000000 ' 



181 



3538944 
67647506447 

25684239974400000 ' 



(27) 



which matches with the result given in Ref. 



28 



The first nontrivial "physical" sector for the TFIM is 
the one-magnon sector, which is a peculiar case of g = 4, 
since the four antiferromagnetic bonds must have a rela- 
tive position similar to the one shown in Fig. 2. We com- 
puted the dispersion £4(^2,, ky) up to order 12. The list of 
all relevant hopping amplitudes is given in Appendix A 1 . 
Noting that e4(fc^, ky) is minimal for (k^, ky) = (0, 0), one 
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FIG. 13: (Color online) Inversion of relations shown in 
Fig. 11. The exponent sbt (for subtracted) refers to the sub- 
cluster subtraction. 



gets the one-magnon gap 



A4 



, 82873487 



19993 , R 

n 

27648 

1901437203257^10 



79626240 1146617856000 
64764934458802909 

23115815976960000 ' *- ' 



Note that our results also match those given in Refs. 19 
and 2c^ up to some trivial rescaling of the Hamiltonian's 
parameters. 

The next "physical" sector for the TFIM is a subspace 
of the 6qp sector which corresponds to the bound states 
discussed above and is illustrated in Fig. 3. In this sector, 
the effective Hamiltonian describes the hopping of the 
bound state, whose center of mass lives on the square lat- 
tice formed by the middles of the bonds (dots in Fig. 4). 
However, one must state that there are two different types 
of sites in this lattice (see Fig. 4) since bonds can be verti- 
cal or horizontal. This is obvious in Fig. 3 where the two 
bound states involve a different pattern of antiferromag- 
netic bonds. Series expansion of the hopping amplitudes 
are given in Appendix A 2 up to order 12. From these 
amplitudes it is clear that the bound state does not have 
the same probability to hop in a given direction, for in- 
stance X, if it is on a horizontal or on a vertical bond. Of 
course, this important distinction also holds for the XXZ 
model discussed in the next section. 

Therefore, one is led to diagonalize a 2 x 2 matrix for 
each value of the center of mass momentum k = (k^, ky). 
One can in particular extract the gap, which is found to 
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FIG. 14: (Golor online) Inversion of relations shown in 
Fig. 12. 



be the lowest of the two energies at fe = (0, 0), and reads 



A- 



275 



11521 



96 27648 
1459322986427 



10 



143327232000 



16400551 , o 

H h^ 

7962624 

101780777359633847 
28894769971200000 



(29) 



,12 



The energy of the second mode at zero center of mass 
momentum is higher and given by 



A+ = 



96 1024 



115 _ 4956689 

39813120 ^ ' 



1720028423 , ^ 880952915946869 ,,2 



17915904000 



9631589990400000 



We do not analyze the associated bound state in detail, 
but just note that this mode decays into the two-magnon 
continuum well before the critical point. Physical impli- 
cations of such a decay will be discussed for the XXZ 
model below. 



B. Exact diagonalization 

To have some finite-size crosschecks of the valid- 
ity of the perturbative expansions and the mapping 
to the Kitaev-type Hamiltonians, performing exact- 
diagonalization (ED) studies on the microscopic models 
defined above demands particular adjustment. First, the 
appearing interactions are of four-body type for TFIM 
in Eq. (10) and of seven-body type for XXZ in Eq. (11). 
Second, the amount of symmetries is enormous, which 
enables us to go to big systems as up to = 50 sites. 
This, however, demands a suitable way to generate the 
symmetrically reduced basis, as a loop through all possi- 
ble spin configurations would be impracticable. To over- 
come the first problem, we apply the Kernel sweeping 
method to efficiently implement many-body interaction 
terms, where details are elaborated on in Ref. 45. Ba- 
sically, a small m-site Hamiltonian of the m-body in- 
teraction is defined, and implemented to sweep over the 
whole lattice basis. To address the second point, we do 
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not loop over all possible TV-site configurations and sub- 
sequently sort out by symmetry constraints, but start 
with an allowed state of the basis and iteratively act on 
it by the respective Hamilton operator. In each step, 
new unprecedented scattering states are collected and 
added to the yet incomplete basis, until a new action 
of H does not produce new scattering states anymore. 
In doing so, as H commutes with the Z2 symmetries 
Eqs. (12) and (14), we make sure that the iteration proce- 
dure only acts within the Z2 charge subsector we want to 
consider. For N = 50, this yields a subblock dimension 
of 2^"/^"^ = 16777216 as alluded to previously, which is 
already in suitable range for Lanczos diagonalization al- 
gorithms. Still, from there, we can further exploit lattice 
symmetries, such as the translational invariance of the 
models, to specify the point of the Brillouin zone we want 
to study. For the Hamiltonians considered, we compute 
the low-energy spectra for clusters up to = 50 sites. As 
shown in Fig. 15, the finite-size corrections do not allow 
us to adequately describe the vicinity of the intermediate 
and large h limit, as the proximity of the continuum and 
the bound-state modes becomes comparable to the finite- 
size splitting scale. Apparently, Fig. 15 explicates that 
the ED data for the largest available system size corre- 
sponds to the expansion data up to /i ^ 0.6. However, in 
the regime where the perturbative expansion data like- 
wise starts to fluctuate a lot depending on the order of 
the expansion, the ED data cannot be used for suitable 
analysis. In the following, it is thus implicitly assumed 
that we use exact diagonalization to crosscheck the imple- 
mentation of the perturbative expansions, but can only 
rely on the latter to highest order to study the interest- 
ing regimes where the bound-state modes approach the 
continuum. 



C. Analysis and gap ratio at the critical point 

As already discussed in Ref. 2N, the series expansion 
for A4 is strongly diverging so that one has to per- 
form some resummation to extract gap values for fi- 
nite h. This also seems to be the case for Ag", al- 
though we only have coefficients up to order 12 for 
this quantity. Unfortunately, standard Pade-type re- 
summation procedure using the order 20 expansion for 
A4 given in Ref. 2S leads to a critical point which 
is still far from the most accurate value. Indeed, 
naive Fade approximants ([10, 10], [8, 12], [12, 8]) lead 
to he = (1.65045, 1.57029, 1.73306), respectively, whereas 
the position of the critical point computed from series ex- 
pansion in the high-field phase"' is he = 1.5216(6). Note 
that for its classical counterpart which is the critical tem- 
perature, higher precision of order 10~^ has been reached 
(see, for instance, Ref. 46 and references therein). 

Here, our aim is to check a result coming from 
field-theoretical calculations performed by Casellc et 
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FIG. 15: (Color online) Thin lines represent the bare ratio 
A(^/A4 for different maximal orders 6, 8, 10, and 12 as a 
function of the magnetic field h. Thick lines correspond to 
different approximations of this ratio. Filled circles denote ED 
data of the mapped model (10) with A*" = 50 sites (bonds of 
the original lattice). Dashed vertical (horizontal) line marks 
the value of the magnetic field (ratio) at the critical point as 
obtained from the numerical diagonalizations of Ref. 17. 



diagonalization"^', giving a ratio of 1.84(3). 

Given that the transition point is defined by A4 = 0, 
a finite ratio at the critical field means that (i) = 
at this point, and (ii) AjT vanishes with the same critical 
exponent v. It can therefore be expected that the direct 
extrapolation of Ag" is at least as complicated as the one 
for A4 which is actually the case. By contrast, one may 
hope that the ratio A^/A4, on which we focus below, 
has a better behavior. 

The bare series of the ratio Ag"/A4 up to order 12 
reads 



^6 
A4 



3 9 ^2 
2 + 16 

156729359 
+ 637009920' 
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768 221184 

27593405457803 
9172942848000 ' 



(31) 
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415396528829457211 
924632639078400000' 



,12 



16,17 



who predict that Ag / Ai\h=ha 



1.8. Re- 



cently, this prediction has been improved using numerical 



and is shown in Fig. 15. Clearly, the bare series is still 
alternating and the convergence is rather poor (h < 0.5) 
so extrapolation schemes are mandatory. First we tried 
standard Fade approximants. We found that all approx- 
imants [n, m] with n + m < 10 give no useful result in the 
sense that the approximant either has a spurious pole or 
shows a diverging behavior well before the critical field 
he- Looking at the Fade approximants with n + m = 12 
we found only two valid cases, namely, [8,4] and [6,6]. 
Still no converging picture emerges because the [8, 4] ap- 
proximant displays a diverging behavior. However, the 
approximant [6, 6] is the first one to behave smoothly 
(see Fig. 15). We conclude that the Fade analysis gives 
no convergent picture and it seems that one needs at least 
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the order 12 scries to catch the physics of this ratio close 
to the critical point. 

To proceed further, we used DlogPade extrapolation 
which is usually more reliable than Pade extrapolations 
for positive quantities. Among all approximants with n > 
4 and m > 4, DlogPade [4, 6] is the only one that has no 
spurious pole. This extrapolation is shown in Fig. 15. It 
can be seen that the ratio seems to behave almost linearly 
as a function of the field close to the critical point. The 
ratio at the critical field he is 1.81 which is very close 
to the numerical value ' . So one finds again that one 
needs at least order 12 to capture the expected behavior. 
Clearly, higher orders are expected to give more valuable 
insights in this quantity but these are beyond the scope 
of this work. 

The relevance of higher orders can be also understood 
from the fact that fluctuations on rather large length 
scales are required to follow this ratio up to the critical 
point. This is in agreement with ED results, displayed 
in Fig. 15, which are qualitatively similar to the bare 
PCUTs series. 



VI. RESULTS FOR THE XXZ MODEL 



A. Series expansion of the low-energy spectrum 



the same steps as for the TFIM discussed in the previous 
section. Series expansions of the hopping amplitudes for 
the bound state are given in Appendix B 2 up to order 
8. One has to diagonalize a 2 x 2 matrix for each value 
of the center-of-mass momentum k = (kx,ky). The gap 
is found at fc — (0,0), and reads 



A7 = 



10 2 323 4 1435321 g 
^"T^ + 540 ^ "^24000"^ 



3809941658983 
' 320060160000 



(34) 



The other bound mode at higher energy with the same 
momentum k = (0, 0) has the expansion 
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(35) 



We would like to point out that the two different orien- 
tations of the bound state are missed in Ref. 15. As a 
consequence, the dispersion for arbitrary momentum of 
the bound state is not correct in this reference. Only the 
gap in this sector matches with our results (which cer- 
tainly means that all hopping amplitudes in Ref. 15 are 
correct). 



Similarly, we computed the low-energy spectrum of 
the XXZ model by the PCUTs method. As previously, 
we restrict the discussion again to the "physical" sub- 
spaces q e {0,4,6} corresponding to the ground state, 
to onc-magnon states and to two-magnon bound states. 
We derive the effective Hamiltonian up to order 10 in 
the OQP sector and up to order 8 in the IQP and 2QP 
sectors. 

For q — 0, we obtain the ground-state energy per dimer 
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(32) 
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which matches with the result given in Ref. 29. 

The onc-magnon sector corresponds again to the first 
nontrivial "physical" sector with g = 4. The four antifer- 
romagnetic bonds must remain in a closed-pack relative 
position such that they share a site in the original square 
lattice (see Fig. 2). The dispersion can be obtained from 
the list of hopping amplitudes given in Appendix B 1. 
The minimum of the dispersion is found at fc = (0,0), 
and one gets for the one-magnon gap, in accordance with 
Ref. 30, 



3 216 
124898889761701 
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(33) 



We now switch to the q = 6 sector that corresponds to 
the two-magnon bound state. Its analysis follows exactly 



B. Fate of the bound states 

An interesting question is to determine when the two- 
magnon bound states decay as a function of A. To 
this end, we first focus on the case of total momen- 
tum k ~ (0,0), describing low-energy physics. Bound 
states decay for this momentum at latest at the Heisen- 
berg point A = 1, where the one-magnon gap closes and 
therefore all multimagnon continua have a gapless spec- 
trum. 

The lowest energy of the two-magnon continuum is 
found at total momentum k — (0, 0) and is given by 
twice the one-magnon gap A4. Bound states start to ac- 
quire a finite lifetime once their energy is degenerate with 
the lower band edge of the two-magnon continuum. Con- 
sequently, we determine the values of A for which ratios 
2A4/A^ are equal to one. In the following, we restrict 
the discussion to the PCUTs results since reliable ED re- 
sults would require enormous system sizes (with 36 sites 
one would only capture terms of order A^). The bound 
state is indeed an extended object and the A term in the 
XXZ Hamiltonian (6) involves nearest-neighbor interac- 
tions, contrary to the purely local h term in the TFIM 
Hamiltonian (1). 

Obviously, the high-energy mode A^ decays first. Us- 
ing different extrapolation schemes such as Pade and 
DlogPade, we find that it disappears for A ~ 0.5401(1), 
i. e., for a rather small value which explains the high ac- 
curacy. Beyond this point, the bound state gains a finite 
lifetime and, strictly speaking, a perturbative derivation 
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FIG. 16: (Color online) Thin blue line represents the bare 
ratio 2A4/A^ for the maximal order 8 as a function of the 
anisotropy parameter A. Thicl: lines correspond to differ- 
ent approximations of this ratio. The dashed horizontal line 
marks the value 1 of the ratio. The value A = 1 corresponds 
to the Heisenberg point. 



of a block-diagonal Hamiltonian becomes impossible"^"'^. 
One can expect that the decayed bound state shows up as 
resonances inside the continua of dynamical correlation 
functions. The two-magnon peak observed in the theo- 
retical Raman response at the Heisenberg point A = 1 
(Refs. 49 and 50) and experimentally detected in the un- 
doped cuprate compounds^'^' '- might be a remnant of 
this bound state. However, a correct physical descrip- 
tion of the decay process is beyond any series expansion 
study. This scenario is further confirmed by the fact that 
any Fade extrapolation of the energy has poles in the 
denominator. 

Concerning the low-energy properties, in analogy to 
the Ising case studied in the previous section, the fate 
of the low-energy mode A^ is more interesting but also 
more challenging. The bare ratio 2A4/A^ together with 
different approximants are shown in Fig. 16. We observe 
that no pole shows up in these approximants. Consis- 
tently, all approximations indicate that the decay takes 
place very close to the Heisenberg point (A ~ 0.97). Tak- 
ing into account that the series have been obtained up 
to order 8 only and have only even orders, one cannot 
tell precisely whether the merging point is located ex- 
actly at the critical point (A = 1), just as in the TFIM, 
or not. One argument in favor of the former scenario 
is that, for any finite-size system, the ground state of 
the SU(2)-invariant Heisenberg point is a singlet and the 
gapped elementary excitation is a threefold degenerate 
triplet, with total magnetization Sz € { — 1,0,-1-1}. The 
excitations with Sz G { — 1,-1-1} can be identified with 
one-magnon excitations found for A < 1, whereas the 
excitation with Sz = has to be a two-magnon bound 
state. The gap of these excitations should furthermore 
vanish in the thermodynamical limit considered in the 
present FCUTs study. 
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FIG. 17: (Color online) Square lattice on which the toric code 
Hamiltonian (36) is defined with periodic boundary condi- 
tions. Sites are represented with dots. We also show a vertex 
s, a plaquette p, a diagonal contour Ci, and an antidiagonal 
contour C2, used to define various operators (see text). 



VII. THE TORIC CODE MODEL IN A 
TRANSVERSE FIELD WITH Jp = 

As mentioned in Sec. II D, the toric code Hamilto- 
nian (10) deserves to be analyzed on its own. Indeed, 
in the original toric code model" ' , Kitaev focused on the 
isotropic coupling Jg = Jp [see Eq. (13) for notations]. 
Here, the bond description of the TFIM leads us to con- 
sider a different situation where (i) Jp = and (ii) a 
magnetic field in the z direction is introduced. Let us 
stress that such a model is very close to the Xu-Moore 
model'" and, in some sense, very similar to the parallel- 
field problem discussed in Refs. 32 and 26,33,34. Most 
importantly, it is exactly the model introduced by Weg- 
ner in his seminal paper '' (see also Ref. ■)■>) but, here. 
Bp operators are conserved quantities that can be in any 
configuration. This crucial difference raises several ques- 
tions that we shall address in the present section. 

Let us first rewrite the Hamiltonian and various oper- 
ators of Sec. II D in the toric code language. The Hamil- 
tonian reads 

iJ = -J^^,-/i^ar, (36) 

s i 

where the spins live on the bonds of a square lattice 
(see Fig. 17) and the As — Jlies '^f operators involve 
the four sites around a vertex s. The plaquette opera- 
tors Bp = Higp '^i ^'"S conserved. For a system defined 
with periodic boundary conditions the cycle operators 
riiec '^i defined on diagonal or antidiagonal contours, 
such as those shown in Fig. 17, are conserved as well. 

The main difference with the Xu-Moore Hamiltonian'^'^ 
is that As operators only act on vertices of the square 
lattice and not on plaquettes. Furthermore, in the Xu- 
Moore model, the cycle operators are still conserved, con- 
trary to the Bp operators. The dynamics of the quasi- 
particles is thus expected to be more constrained than in 
the Xu-Moore model which already exhibits dimensional 
reduction. 
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In the following, we shall discuss separately the two 
limits (large and small J/ h) for which we computed per- 
turbatively the low-energy spectrum. One can already 
note that in the small- field limit, the system is in a topo- 
logical phase (i. e., the ground-state degeneracy depends 
on the surface genus) whereas at large field the ground 
state is obviously unique. Thus, one expects a (topolog- 
ical) quantum phase transition when varying the ratio 
J/h. 



A. Large- field limit J 

For J = 0, elementary excitations are usual magnons 
obtained by flipping any spins from the ground state 
which is the state fully polarized in the field (z) direc- 
tion. PCUTs formalism will then give a dressed magnon 
description when switching on J-''. Setting h = 1/2, the 
ground-state energy per bond (which coincides with the 
Oqp level) is obtained from Eq. (27) by replacing h by J. 

However, for this model, arbitrary qp sectors are al- 
lowed although conservation of -Bp's imposes severe con- 
straints on the spectrum. For q = I, one only has a 
nondispersive (excited) level at energy 



. J2 3J4 

Ai = 1 \ 

2 32 

2014178639J 



31J^ 



299233J** 



768 15925248 



10 



764411904000 ' ^^'^^ 
since any displacement of this localized dressed magnon 
would modify the Bp's configuration. 

For q = 2, the only possible dynamics is provided by 
some flip-flop processes depicted in Fig. 18 (we do not 
consider the case where the two dressed magnons are far 
apart). Both configurations shown in Fig. 18 yield the 
same gap. which reads 

13J4 



Ao = 2 



29J=' 



445 



192 768 13824 
608839J*^ 2462069J9 



739J^ 
' 36864 ~ 79626240 
21097903J1" 



9555148800 



' 152882380800' 
Note that in the Xu-Moore model, the pair of nearest- 
neighbor magnons (top of Fig. 18) is allowed to hop at 
any distance (but in a one-dimensional stripe) 

We have not computed energies of states with more 
magnons. Let us simply mention the following two 
points. For q = 3, the problem is similar to g = 2 and 
no dynamics of the QPs is allowed apart from flip-flip 
processes. The first dispersive mode is the one involving 
four magnons, just as in the Xu-Moore model^ ^ 



B. Small-field limit < J 

Let us now turn to the opposite limit where J is much 
larger than h and for which the system is in a topological 
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FIG. 18: (Color online) Illustration of the two-magnon dy- 
namics in the high- field limit of Hamiltonian (36). Dressed 
magnons are depicted as crosses. The + and — signs refer to 
the eigenvalue -1-1 or —1 of the plaquette operators Bp. 



phase. For convenience, let us set J = 1/2. Since the Bp 
operators do not appear in the Hamiltonian, all states 
are macroscopically degenerate when h = 0. It is thus 
natural to wonder what the effect of the perturbation is. 
A calculation at fourth order of the effective Hamiltonian 
in the sector where there is no effective vertex excitation 
(still using the PCUTs formalism) yields 



H^ff = - + — + — -] ^> Br, 

1428/ 2^^ 



(39) 



where N is the number of sites (and thus twice the num- 
ber of plaquettes or of vertices). We thus see that, once 
the magnetic field is switched on, the ground state be- 
longs to the sector where each operator Bp has eigenvalue 
+1. As a check, note that setting all Bp operators to 
one in the above formula gives the ground-state energy, 
which matches with Eq. (8) of Ref. '2(> (apart from a triv- 
ial constant shift since no Bp operator appears in our 
Hamiltonian, and after setting to zero in this equa- 
tion). The fact that the sector with no plaquette exci- 
tation is selected by the perturbation is consistent with 
an Ising-type quantum phase transition. Indeed, in this 
sector, the model with Hamiltonian (36) is dual to the 
two-dimensional transverse-field Ising model (which fol- 
lows from the same arguments as those used in Ref. '-yi). 



VIII. CONCLUSIONS 

We have studied two-magnon bound states in the 
TFIM and XXZ models using PCUTs around the Ising 
limit. This method gives a nice physical picture of these 
excitations in terms of dressed magnons which can be 
of two kinds on a square lattice depending on the "frus- 
trated bond" localizing it (horizontal or vertical) . Conse- 
quently, one obtains two different bound modes contrary 
to previous claims'^ For both systems, we find that 
the lowest-energy gap vanishes at the critical point, or 
at least in a neighborhood that cannot be distinguished 
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from this point with actual series expansions. Solving 
this issue, especially for the XXZ model, would clearly 
require much higher orders or a nonperturbative treat- 
ment as the one considered in Refs. IG and 17 for the 
TFIM. 

From a methodological point of view, we have adapted 
Entings' finite-lattice method commonly used in statis- 
tical mechanics to quantum problems. This approach 
basically consists in a cluster embedding which consider- 
ably increases the efficiency (from the time and memory 
point of view) of the PCUTs method we used here. Note 
that for the XXZ model, such an improvement allowed 
us to reach the same maximum order as standard se- 
ries expansions techniques (see Refs. 19 and 15) based 
on a more sophisticated graph analysis. Clearly, adapt- 
ing such a graph description to PCUTs is a crucial issue 
which is currently under study and should allow one to 
reach higher orders. 
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Appendix A: Low-energy spectrum of the TFIM 

1. One-magnon hopping amplitudes 

One-magnon excitations are located on sites of the 
square lattice. The corresponding hopping amplitudes 
ti j of i sites in x direction and j sites in y direction are 
given below. Hopping amplitudes that are not given can 
be deduced from the symmetries of the lattice and the 
hermiticity of the Hamiltonian. 
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2. Two-magnon bound state hopping amplitudes 
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The two-magnon bound states are located at the cen- 
ters of the bonds on the square lattice and exist in two 
kinds : the bound state can indeed live on a horizontal or 
a vertical link denoted h and v respectively. The hopping 
amplitudes of a horizontally oriented bound state t^^'" 
are listed below. The corresponding hopping elements 
for the vertically oriented bound state can be deduced 
by reversing the x and y components, i. e., tj^" = t^'"- 
Again, the symmetries of the lattice and the hermiticity 
of the Hamiltonian have been used to restrict the number 
of given amplitudes. 
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Appendix B: Low-energy spectrum of the XXZ 

Notations are the same as in Appendix A. 

1. One-magnon hopping amplitudes 
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2. Two-magnon bound state hopping amplitudes 
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